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Abstract 

The southeastern US has experienced recurrent drought during recent decades. Increasing demand for water, as 
precipitation decreases, exacerbates stress on the aquatic biota of the Southeast: a global hotspot for freshwater mussel, 
crayfish, and fish diversity. Freshwater unionid mussels are ideal candidates to study linkages between ecophysiological and 
behavioral responses to drought. Previous work on co-occurring mussel species suggests a coupling of physiology and 
behavior along a gradient ranging from intolerant species such as Pyganodon grandis (giant floater) that track receding 
waters and rarely burrow in the substrates to tolerant species such as Uniomerus tetralasmus (pondhorn) that rarely track 
receding waters, but readily burrow into the drying sediments. We utilized a next-generation sequencing-based RNA-Seq 
approach to examine heat/desiccation-induced transcriptomic profiles of these two species in order to identify linkages 
between patterns of gene expression, physiology and behavior. Sequencing produced over 425 million 100 bp reads. Using 
the de novo assembly package Trinity, we assembled the short reads into 321,250 contigs from giant floater (average length 
835 bp) and 385,735 contigs from pondhorn (average length 929 bp). BLAST-based annotation and gene expression 
analysis revealed 2,832 differentially expressed genes in giant floater and 2,758 differentially expressed genes in pondhorn. 
Trancriptomic responses included changes in molecular chaperones, oxidative stress profiles, cell cycling, energy 
metabolism, immunity, and cytoskeletal rearrangements. Comparative analyses between species indicated significantly 
higher induction of molecular chaperones and cytoskeletal elements in the intolerant P. grandis as well as important 
differences in genes regulating apoptosis and immunity. 
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Introduction 

As drought conditions recur with increasing frequency and 
severity in the Southeastern U.S., sessile aquatic organisms in 
freshwater ecosystems, particularly in streams, are bearing the 
brunt of these environmental perturbances [1-3]. Freshwater 
unionid mussel populations are already among the most endan- 
gered groups of organisms in the world [4]. Growing efforts to 
document the stunning diversity and varied life history strategies of 
unionids are also recording the measurable impacts of altered 
stream flows and thermal profiles on survival, recruitment, 
reproductive strategies, and community structure of these species 
[2,3,5-7]. 

The ability of individual species to tolerate drought conditions 
depends on several factors including severity and duration of the 
disturbance, stream habitat (debris, pools, etc.), and differences in 
behavioral and physiological responses to emersion/desiccation 
and heat stress [3]. Indeed, a recent study utilizing field and 



laboratory experiments revealed links between behavioral respons- 
es, physiological tolerances, and survival in three co-existing 
mussel species (Uniomerus tetralasmus (pondhorn); Pyganodon grandis 
(giant floater); and Lampsilis straminea (fatmucket)). There, the 
authors observed that a burrowing response in pondhorn was 
correlated with a higher tolerance to desiccation and higher 
survival (77%) compared to the water-tracking behavior and low 
tolerance and survival (0%) of giant floater [3]. 

The identification of underlying genetic mechanisms regulating 
these behavioral and physiological differences would provide key 
insights into adaptive responses of freshwater mussels to heat stress 
and drought. Previous studies in marine shellfish species have 
explored connections between variations in gene and protein 
expression and differences in latitudinal adaptation, differential 
success of native and invasive species, and resistance to summer 
mortality in the context of heat stress [8-13]. While such studies 
have traditionally required time-consuming and expensive gener- 
ation of molecular resources and have, therefore, been limited to a 
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handful of key species, new RNA-Seq approaches have dramat- 
ically expanded our ability to carry out transcriptome profiling in 
non-model species [14]. Indeed, using RNA-Seq we recendy 
identified the components of a classical heat stress response in the 
unionid Villosa lienosa [15]. Of even greater interest than highly 
conserved stress factors and responses [16], however, may be the 
identification of divergent responses and gene components 
accounting, at least in part, for heightened drought tolerance in 
some species or populations. Recent sequencing of the Pacific 
oyster genome revealed a starding diversity of genes available for 
adaptation to environmental stress, including 88 heat shock 
protein 70 (HSP70) genes (compared to 17 in humans) and 48 
inhibitor of apoptosis (LAP) genes (compared to 8 in humans). The 
finding not only highlights the importance of these gene families in 
molluscs but also emphasizes the need for global expression 
profiling of particular species of interest rather than examination 
of a conserved, narrow subset of heat shock genes based on 
mammalian paradigms. Here, we conducted RNA-Seq-based 
transcriptome profiling on P. grandis and U. tetralasmus exposed to 
experimental heat stress/ desiccation. Our main objectives were: (i) 
to compare the ability of Trinity and Trans-ABySS to assemble 
high quality transcriptomes for both species; (ii) identify key shared 
and divergent responses to drought in unionid mussels; and (iii) 
determine whether there were consistent differences between the 
tolerant and sensitive unionid species in both individual gene and 
pathway level responses to drought stress. The information gained 
here may serve as a foundation to guide natural resource managers 
in assessing the adaptive potential of freshwater mussel species in 
the face of climate change and should aid in the development of 
molecular tools to monitor stress levels of mussels in drought- 
stricken streams. 

Materials and Methods 

Experimental animal and tissue sampling 

Two species of freshwater mussel, P. grandis and U. tetralasmus, 
were utilized in experiments. Pyganodon grandis have been desig- 
nated a species of lowest conservation concern, and U. tetralasmus a 
species of moderate conservation concern in Alabama. Uniomerus 
tetralasmus were collected from Opintlocco Creek, located in the 
Tallapoosa catchment of east central Alabama, with permission 
from the J. W. Huskey family and under scientific collection 
permits from the Alabama Department of Conservation and 
Natural Resources. Pyganodon grandis were produced in experimen- 
tal ponds at the South Auburn Fisheries Research Station, Auburn 
University as part of previous experiments by the Crustacean and 
Molluscan Ecology Lab. Prior to experiments, all mussels were 
submerged in a large water bath and maintained at 21°C. Pond 
water was used throughout the study so mussels could feed on 
natural food sources. All mussels were placed individually in 
500 mL plastic cups within the water bath. Each cup was filled 
with sand and mussels were allowed to position themselves 
naturally in the substrate. After three days of acclimation at 21°C, 
temperature was increased 3°C/d over 4 d to 33°C for the 
experimental mussels, while temperature (21°C) and water level 
remained constant for the duration of the experiment for the 
control mussels. Upon reaching 33°C, mussels were maintained at 
that temperature for 2 d. On day 3 at 33°C, the water level in the 
water bath was lowered below the top of the cups and the water in 
each cup was drained. 

At 24, 48, and 72 h after initiation of desiccation/emersion, 
non-lethal foot tissue samples were collected from 15 individuals 
per species using biopsy forceps. Mussels were randomly assigned 
to sampling time groups. After the 72 h samples were collected, 



the recovery phase of the experiment began. The remaining 
experimental mussels were submerged and temperatures were 
dropped to 21°C by 3°C/d for 4 d. After one day at 21°C, 15 
mussels per species were sampled (5 d recovery timepoint). Tissues 
were flash frozen with liquid nitrogen immediately upon collection 
and stored in a — 80°C freezer until RNA extraction. 

RNA extraction, library construction and sequencing 

From each timepoint/treatment group, three replicate pools of 
5 individual samples/pool were constructed. Pooled tissue samples 
were ground to a fine powder in the presence of liquid nitrogen. 
RNA was extracted from ground tissue using an RNeasy Kit 
(Qiagen, Valencia, California) following the manufacturer's 
instruction. RNA concentration and integrity of each sample were 
measured on an RNA 6000 Nano Bioanalysis chip using the 
Agilent Bioanalyzer 2100 on the RNA 6000 Nano Bioanalysis 
chip, only the samples with R.I.N, value >8 were used for the 
sequencing. 

The control group (submerged 21°C) and 72 h (desiccation 
33°C) experimental groups from each species were selected for 
RNA-Seq library construction, Illumina sequencing and expres- 
sion analysis. 

RNA-Seq library preparation and sequencing were carried out 
by HudsonAlpha Genomic Services Lab (Huntsville, AL, USA). 
cDNA libraries were constructed with 2.14-3.25 ug of starting 
total RNA and utilized the IlluminaTruSeq RNA Sample 
Preparation Kit (Illumina), as detailed in the TruSeq protocol. 
The libraries were amplified with 15 cycles of PGR and contained 
TruSeq indexes within the adaptors, specifically indexes 1-12. 
Finally, amplified library yields were 30 ul of 19.8-2 1 .4 ng/ ul with 
an average length of — 270 bp, indicating a concentration of 1 10- 
140 nM. After KAPA quantitation and dilution, the libraries were 
clustered 12 per lane and sequenced on an Illumina HiSeq 2000 
instrument with 100 bp PE read chemistry. 

De novo assembly of sequencing reads 

The de novo assembly of short reads was performed on the mussel 
short reads using both ABySS and Trinity [17,18], versions 1.3.2 
and the 2012-10-05 editions, respectively. Before assembly, raw 
reads were trimmed by removing adaptor sequences and 
ambiguous nucleotides. Reads with quality scores less than 20 
and length below 30 bp were all trimmed. The resulting high- 
quality sequences were used in the subsequent assembly. 

In ABySS, briefly, the clean reads were first hashed according to 
a predefined k-mer length, the 'k-mers'. After capturing overlaps 
of length k-1 between these k-mers, the short reads were assembled 
into contigs. The k-mer size was set from 50 to 96, assemblies from 
all k-mers were merged into one assembly by Trans-ABySS 
(version 1.4.4). 

In Trinity, briefly, the raw reads were assembled into the unique 
sequences of transcripts in Inchworm via greedy k-mer extension 
(k-mer 25). After mapping of reads to Inchworm contigs, Chrysalis 
incorporated reads into de Bruijn graphs and the Butterfly module 
processed the individual graphs to generate full-length transcripts. 

In order to reduce redundancy, the assembly results from 
different assemblers were passed to CD-Hit version 4.5.4 [19] and 
CAP3 [20] for multiple alignments and consensus building. The 
threshold was set as identity equal to 1 in CD-Hit, the minimal 
overlap length and identity equal to 100 bp and 99% in CAP3. 

Gene annotation and ontology 

All contigs in the resulting Trinity assembly were used as queries 
for BLASTX searches against the UniProtKB/SwissProt database 
and the National Center for Biotechnology Information (NCBI) 
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non-redundant (nr) protein database. The cutoff Expect value (E- 
value) was set at le— 5 and only the top hit result was assigned as 
the annotation for each contig. Gene ontology (GO) annotation 
analysis was processed using the UniProtKB/SwissProt results in 
Blast2GO (version 2.5.0), which is an automated tool for the 
assignment of gene ontology terms. The final annotation file was 
generated after gene-ID mapping, GO term assignment, annota- 
tion augmentation and generic GO-Slim processes and catego- 
rized with regard to Biological Process, Cellular Component and 
Molecular Function. 

Identification of differentially expressed contigs 

Differentially expressed contig analysis was performed using the 
RNA-Seq module and the expression analysis module in CLC 
Genomics Workbench. The high quality reads from each sample 
were mapped onto the Trinity reference assembly using CLC 
Genomics Workbench software. At least 95% of the bases were 
required to align to the reference and no more than two 
mismatches were allowed. The total mapped reads number for 
each transcript was determined and then normalized to detect 
RPKM (Reads Per Kilobase of exon model per Million mapped 
reads). After scaling normalization of the RPKM values, fold 
changes were calculated [21]. The proportions-based test was used 
to identify differentially expressed contigs between the control 
group and heat/desiccation group with three replications in each 
group with the threshold for significant differences set at a 
corrected p-value of /><0.05 [22]. Transcripts with absolute fold 
change values greater than 1.5 were defined as differentially 
expressed genes, which were divided into up- and down- regulated 
groups for further analysis. 

Differentially expressed contigs were also used as queries for 
BLASTX searches against the Crassostrea gigas (Pacific oyster) 
protein database from NCBI with the cutoff Expect value (E-value) 
of le— 5. Contigs from P. grandis and U. tetralasmus with shared 
annotation (nr and/ or C. gigas) based on BLAST analysis were 
subjected to additional analyses to establish orthology including 
reciprocal BLAST, sequence alignment, and phylogenetic analysis. 
Differentially expressed contigs from a given species with no clear 
orthologous gene(s) in the other species were regarded here as 
unique. 

Gene ontology and enrichment analysis 

GO analysis and enrichment analysis between significandy 
expressed GO terms and overall reference assembly GO terms was 
processed to select overrepresented GO annotations in the 
differentially expressed genes using Ontologizer 2.0 using the 
Parent-Child-Intersection method with a Benjamini-Hochberg 
multiple testing correction [23,24]. GO terms for each gene were 
obtained by utilizing UniProtKB/SwissProt database annotations 
for the unigene set. The difference in the frequency of GO terms 
annotation in the differentially expressed genes sets were 
compared to the overall mussel reference assembly for each 
species. The threshold was set as FDR value <0.1. 

Functional groups and pathways encompassing the differentially 
expressed genes were identified based on GO analysis, pathway 
analysis based on the Kyoto Encyclopedia of Genes and Genomes 
(KEGG) database, and manual literature review. Based on these 
analyses, key genes were organized into broad functional 
categories. Fold change differences between species for a given 
gene in each category were assessed for significance using a simple 
t-test on log-transformed RPKM values. Category-level differences 
in fold changes between species were assessed using a paired t-test. 



Experimental validation — QPCR 

Sixteen significantly expressed genes, including 7 shared genes, 
4 unique genes of P. grandis and 5 unique genes of U. tetralasmus, 
with different expression patterns were selected for validation of 
the RNA-Seq results using real time QPCR. Beta-actin was set as 
the reference gene in both species. Primers were designed based on 
contig sequences using Primer Premier 6 (Table SI). Tissue 
samples from the 72 h experimental group and the control group 
were used as the template for QPCR. cDNA was synthesized using 
the iScript cDNA Synthesis Kit (Bio-Rad) according to manufac- 
turer's protocol. The iScript chemistry uses a blend of oligo-dT 
and random hexamer primers. All the cDNA products were 
diluted to 250 ng/u.1 and utilized for the quantitative real-time 
PCR reactions using the SsoFast EvaGreen Supermix on a CFX96 
real-time PCR Detection System (Bio-Rad Laboratories, Hercules, 
CA). The thermal cycling profile consisted of an initial denatur- 
ation at 95°C (for 30 s), followed by 40 cycles of denaturation at 
94°C (5 s), and an appropriate annealing/extension temperature 
(58°C, 5 s). An additional temperature ramping step was utilized 
to produce melting curves of the reaction from 65 °C to 95 °C. 
Results were expressed relative to the expression levels of beta 
actin in each sample using the Relative Expression Software Tool 
(REST) version 2009 [25]. The biological replicate fluorescence 
intensities of the control and treatment products for each gene, as 
measured by crossing-point (Ct) values, were compared and 
converted to fold differences by the relative quantification method. 
Expression differences between groups were assessed for statistical 
significance using a randomization test in the REST software. The 
mRNA expression levels of all samples were normalized to the 
levels of beta actin gene in the same samples. Test amplifications 
were conducted with ensure that beta actin and target genes were 
within an acceptable range. A no-template control was run on all 
plates. QPCR analysis was repeated in triplicate runs (technical 
replicates) to confirm expression patterns. 

Temporal expression analysis of key genes 

A subset of key classical heat shock proteins and species-specific 
genes differentially expressed at 72 h were selected to examine 
their broader transcriptional profiles. Primers, reference gene, and 
reaction conditions were the same as described above. Relative 
expression profiles were captured for the experimental samples of 
each species at the 24 h, 48 h, and 5 d recovery timepoints. Data 
were analyzed as described above. 

Results 

Sequencing of short expressed reads from mussel foot 
tissue 

Illumina-based RNA-sequencing (RNA-Seq) was carried out on 
RNA extracted from replicated, pooled foot tissue samples from P. 
grandis and U. tetralasmus exposed to control and heat stress/ 
desiccation (72 h) experimental conditions. Reads from different 
samples were distinguished through the use of multiple identifier 
(MID) tags. A total of 427.5 million reads were generated on an 
Illumina HiSeq2000 instrument, including 201.1 million reads 
from P. grandis, and 226.4 million reads from U. tetralasmus (Table 
S2). At least 29 million reads were generated from each barcoded 
sample with an average of 35.6 million reads/library. Raw data 
was archived at the NCBI Sequence Read Archive (SRA) under 
Accession SRP026193 . 

De novo assembly of mussel transcriptomes 

De novo assembly of RNA-Seq reads can be achieved using 
several assembly algorithms and software programs. We had 
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previously developed an in-house bioinformatics pipeline around 
Trans-ABySS and demonstrated superior performance relative to 
other assembly options in several species including freshwater 
mussel [15,26-28]. However, the recently developed Trinity 
software package provides several potential advantages for 
transcriptome profiling in non-model species [18]. Therefore, 
here we sought to compare critical metrics of performance in 
transcriptome assembly between Tran-ABySS and Trinity. 

Trans-ABySS. From P. grandis cleaned reads, Trans-ABySS 
generated 804,906 contigs including 178,956 contigs longer than 
1000 bp. Average contig length was 772.6 bp while N50 was 
801 bp. After filtration of redundant contigs by CD-Hit and 
CAP3, 56.6% of total contigs remained (455,659) with average 
length of 616.0 bp. The percentages of reads mapped in pairs and 
mapped to the final reference were 62.8% and 83.7%, respectively 
(Table 1). 

From U. tetralasmus cleaned reads, Trans-ABySS generated 
499,413 contigs including 147,046 contigs longer than 1000 bp. 
Average contig length was 963.0 bp while N50 was 1,503 bp. 
After filtration of redundant contigs by CD-Hit and CAP3, 52.4% 
of total contigs remained (261,489) with average length of 
885.5 bp. The percentages of reads mapped in pairs and mapped 
to the final reference were 65.9% and 82.0%, respectively 
(Table 1). 

Trinity. From P. grandis cleaned reads, Trinity generated 
336,799 contigs including 66,959 contigs longer than 1,000 bp. 
Average contig length was 970.1 bp while N50 was 2,100 bp. 
After filtration of redundant contigs by CD-Hit and CAP3, 95.4% 
of total contigs remained (321,250) with average length of 
835.0 bp. The percentages of reads mapped in pairs and mapped 
to the final reference were 79.2% and 84.5%, respectively 
(Table 1). 

From U. tetralasmus cleaned reads, Trinity generated 405,996 
contigs including 81,095 contigs longer than 1,000 bp. Average 
contig length was 968.1 bp while N50 was 2,233 bp. After 
nitration of redundant contigs by CD-Hit and CAP3, 95.0% of 
total contigs remained (385,735) with average length of 929.3 bp. 
The percentages of reads mapped in pairs and mapped to the final 
reference were 79.5% and 83.8%, respectively (Table 1). 

Best assembly selection. Comparison of the assemblies 
generated by Tran-ABySS and Trinity (Table 1) made clear that, 
although Tran-ABySS generated a larger initial count of contigs, 
redundancy was much higher than observed with Trinity. CD- 
HIT/CAP3 trimmed over 40% of Trans-ABySS contigs com- 



pared with less than 5% of Trinity contigs. Additionally, Trinity 
produced contigs with greater N50, average length and mapped 
reads metrics. Trinity was able to assemble transcript pieces into 
longer contigs with average lengths over 200 bp greater than 
Trans-ABySS and a high rate of paired read mapping. Consid- 
ering these results, we utilized the Trinity assembly for subsequent 
analysis. Trinity contig assemblies for P. grandis and U. tetralasmus 
are available upon request. 

Gene identification and annotation 

Trinity contigs were used to query the UniProtKFS/SwissProt 
and NCBI non-redundant (nr) protein databases by BLASTX 
(Table 2). 

In P. grandis, 57,469 contigs had significant BLASTX hits 
against 22,616 unique nr proteins. Of these, 6,120 contigs had top 
BLAST hits against hypothetical gene matches (predicted proteins 
and/or those with no annotation). Additionally, 11,877 unigenes 
could be identified based on more stringent criteria of a BLAST 
score SI 00 and E-value S le — 20 (quality matches). The same 
BLAST criteria were used to query the UniProt database. 



Table 2. Summary of gene identification and annotation of 



assembled mussel contigs based on BLAST homology 
searches against various protein databases (UniProt, nr). 






P. grandis 




U. tetralasmus 




UniProt 


NR 


UniProt 


NR 


Contigs with putative gene 
matches 


48,448 


57,469 


56,170 


66,381 


Percentage of annotated contigs 


14.38% 


1 7.06% 


1 3.84% 


16.35% 


Annotated contigs >500 bp 


28,205 


32,290 


33,595 


38,375 


Annotated contigs>1000 bp 


11,773 


13,243 


14,391 


16,017 


Unigene matches 


1 5,045 


22,616 


15,170 


20,610 


Hypothetical gene matches 


0 


6,120 


0 


6,101 


Quality Unigenematches 


11,278 


11,877 


11,324 


10,699 



Putative gene matches were at E-value < le — 5. Hypothetical gene matches 
denote those BLAST hits with uninformative annotation. Quality unigene hits 
denote more stringent parameters, including score>100, E-value < le— 20. 
doi:1 0.1 371 /joumal.pone.0089481 .t002 



Table 1 . Summary of de novo assembly results of I 
ABySS and Trinity. 



lumina cleaned reads from P. grandis and U. tetralasmus foot tissue using Trans- 





P. grandis 




U. tetralasmus 






Trans-ABySS 


Trinity 


Trans-ABySS 


Trinity 


Contigs 


804,906 


336,799 


499,413 


405,996 


Large contigs (>1000 bp) 


178,956 


66,959 


147,046 


81,095 


N50 (bp) 


801 


2,100 


1,503 


2,233 


Average contig length 


772.6 


970.1 


963.0 


968.1 


Contigs (After CD-HIT-EST+ CAP3) 


455,659 


321,250 


261,489 


385,735 


Percentage contigs kept after redundant removing 


56.6% 


95.4% 


52.4% 


95.0% 


Average length (bp) (After CD-HIT-EST+ CAP3) 


616.0 


835.0 


885.5 


929.3 


Reads mapped in pairs (%) 


62.8 


79.2 


65.9 


79.5 


Reads mapped to final reference (%) 


83.7 


84.5 


82.0 


83.8 



doi:1 0.1 371 /joumal.pone.0089481 .t001 
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In U. tetralasmus, 66,381 contigs had significant BLASTX hits 
against 20,610 unique nr proteins. Of these, 6,101 contigs had top 
BLAST hits against hypothetical gene matches (predicted proteins 
and/or those with no annotation). Additionally, 10,699 unigenes 
could be identified based on more stringent criteria of a BLAST 
score &100 and E-value S le — 20 (quality matches). The same 
BLAST criteria were used to query the UniProt database. 

Identification and analysis of differentially expressed 
genes 

Differential expression analyses in comparison to control 
samples were carried out for both the P. grtmdis and U. tetralasmus 
experimental samples (Table S3; Table 3). In P. grandis, a total of 
2,559 genes (unique annotated contigs with significant nr database 
BLASTX identities) including 1 ,499 up-regulated genes and 1 ,060 
down-regulated genes were differentially expressed (S 1.5-fold, 
corrected p-value <0.05) following 72 h of exposure to heat and 
desiccation stress. Similarly, in U. tetralasmus, 2,532 genes 
including 1,594 up-regulated genes and 938 down-regulated genes 
were differentially expressed. Read coverage, critical in accurate 
determination of fold change, averaged 797 and 567 reads/contig 
for P. grandis and U. tetralasmus, respectively. We also used the 
differentially expressed contigs of both species as queries against 
annotated proteins from the recently sequenced Pacific oyster, C. 
gigas genome (the only publicly available genome among bivalve 
species [29]). Annotating based on hits to the C. gigas protein 
database identified 1,863 differentially expressed genes in P. grandis 
and 1,744 differentially expressed genes in U. tetralasmus. 

While fewer contigs were annotated by C. gigas, its larger, 
comprehensive protein database aided in initial identification of 
"shared" genes present in sets of differentially expressed genes 
from both mussel species. For example, we identified .318 and 380 
significantly up-regulated contigs with the same top hit in both 
species using nr and C. gigas annotation, respectively. Similarly, 
190 contigs shared the same top hit and were down-regulated in 
both species using C. gigas annotation. Finally, there are 130 and 
162 differentially expressed genes in each species shared the same 
BLAST identity but showed fold changes in opposite directions. 

Enrichment and pathway analysis 

Gene ontology (GO) annotations were assigned using Blas- 
t2GO. In P. grandis, 7,571 GO terms were assigned to 2,559 
unique genes including 2,660 (35.1%) cellular component terms, 
1,629 (21.5%) molecular functions terms and 3,282 (43.4%) 
biological process terms. In U. tetralasmus, a total of 6,764 GO 
terms were assigned to 2,532 unique gene matches including 2,478 

Table 3. Statistics of differentially expressed genes of 
P. grandis and U. tetralasmus annotated by nr and C. gigas 
protein database following drought challenge. 





P. grandis 




U. tetralasmus 




NR 


C. gigas 


NR 


C. gigas 


Up-regulated 


1,499 


1,093 


1,594 


1,126 


Down-regulated 


1,060 


770 


938 


648 


Total unigenes 


2,559 


1,863 


2,532 


1,774 


Reads per contig 


797 


718 


569 


586 



Values indicate contigs/genes passing cutoff values of fold change 5:1.5 
(corrected p<0.05) and read number per contig >5. Reads/contig refers to 
average contig size. 
doi:1 0.1 371 /joumal.pone.0089481 .t003 



(36.6%) cellular component terms, 1530 (22.6%) molecular 
functions terms and 2,756 (40.8%) biological process terms. The 
percentages of annotated sequences of the two species assigned to 
GO terms are shown in Figure SI. 

The differendy expressed unique genes were then used as inputs 
to perform enrichment analysis using Ontologizer. Parent-child 
GO term enrichment analysis was performed to detect signifi- 
candy overrepresented GO terms. A total of 23 and 50 
overrepresented terms with p (FDR-corrected) <0. 1 were detected 
in P. grandis and U. tetralasmus, respectively. Potentially informative 
higher level GO terms were carried for further pathway analysis. 
These included indications of changes in protein folding in both 
species, and species-specific enrichment of cellular transport, cell 
structure and energy metabolism in P. grandis and genes negatively 
regulating cell death and the immune response in U. tetralasmus 
(Table S4). However, the general lack of functionally annotated 
genes in molluscan species limited the extent of GO analysis that 
could be carried out. 

Based on enrichment analysis, manual annotation, and litera- 
ture searches, representative key genes were organized into 7 
functional categories encompassing important shared and poten- 
tially species-specific response signatures to heat stress/desiccation 
(Table 4). These included chaperone/heat shock proteins, 
antioxidant/ oxidative stress response, cell proliferation/apoptosis, 
energy metabolism, immune, cytoskeletal and uncategorized 
(other). While imputed putative functional roles of these genes 
are covered in-depth below (Discussion), several general trends 
were apparent. First, while a robust heat-shock response was 
mounted in both species, the drought-susceptible P. grandis 
manifested significantly higher and broader up-regulation of 
molecular chaperones at 72 h (p = 0.0277, paired t-test). Second, 
significandy different category-level fold changes were observed in 
genes regulating cytoskeletal arrangements (p = 0.0149), with 
higher up-regulation of cell fiber elements observed for P. grandis. 
Third, individual key genes within each category had significandy 
different responses to drought stress, including, for example, those 
regulating inhibition of apoptosis and immunity (Table 4). 

Validation of RNA-Seq profiles by QPCR 

In order to validate the differentially expressed genes identified 
by RNA-Seq, we selected 16 genes for QPCR confirmation from 
P. grandis and U. tetralasmus (Table S5). Selection of genes was based 
on putative functions in response to heat stress/desiccation and 
pathway analyses. Melting curve analysis revealed a single product 
for all tested genes. Fold changes from QPCR were compared 
with the RNA-Seq expression analysis results. QPCR results for P. 
grandis were significantly correlated with RNA-Seq results (avg. 
correlation coefficient, R = 0.87; p-value <0.001; Figure 1A). 
Similarly, with the exception of calpain 5, QPCR results for U. 
tetralasmus closely matched those observed by RNA-Seq results 
(avg. correlation coefficient, R= 0.89; p-value <0.001; Figure IB). 
Overall, the QPCR results indicated the reliability and accuracy of 
the Trinity reference assembly and the RNA-Seq-based tran- 
scriptome expression analysis. 

Temporal expression of key heat/desiccation signatures 

Given our focus on the 72 h timepoint for RNA-Seq expression 
analysis, we sought to extend the profiles of several key genes to 
earlier (0 h, 24 h, 48 h) timepoints as well as the recovery phase 
(5 d). While financial constraints limited our ability to examine 
these timepoints comprehensively by RNA-seq, we wished to 
examine whether key genes identified at 72 h could serve as 
sensitive markers for stress at other timepoints. We selected four 
heat shock genes that were up-regulated to differing extents in 
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Table 4. Key differentially expressed genes following drought challenge in P. grandis and U. tetralasmus. 



Gene description 


NCBI Accession 
number 


P. grandis 






U. tetralasmus 






t-test 
p<0.05 






Contig ID 


FC 


p-value 


Contig ID 


FC 


p-value 




Chaperone/HSP (p= 0.0277) 


60 kDa heat shock protein 


EKC31 862.1 


fl_ctg_778 


2.7 


4.7E-05 


ph_ctg_636 


2.0 


7.1E-03 




78 kDa glucose-regulated protein 


EKC33663.1 


fl_ctg_1364 


6.0 


3.6E-43 


ph_ctg_1443 


3.8 


2.6E-31 




Alpha-crystallin B chain 


EKC27067.1 


fl_ctg_1 125 


35.8 


2.8E-1 1 


ph_ctg_1374 


7.2 


1.5E-152 


* 


DnaJ-like protein subfamily A member 1 


EKC32228.1 


fl_ctg_1 29 


3.4 


1.3E-02 


ph_ctg_1424 


2.2 


1 .9E-04 




DnaJ-like protein subfamily A member 2 


EKC 18925.1 


fl_ctg_188 


2.0 


2.5E-15 


ph_ctg_58 


1.6 


2.2E-05 


* 


Endoplasmin 


EKC38233.1 


fl_ctg_1 292 


3.1 


7.1E-35 


ph_ctg_790 


2.1 


1.8E-08 




Heat shock 10kDa protein 


NP_001 084708.1 


fl_ctg_489 


5.3 


4.9E-03 


ph_ctg_1163 


2.1 


4.7E-06 


* 


Heat shock 70 kDa protein 1 2A, type 1 


EKC29170.1 


fl_ctg_1046 


-2.7 


1 .5E-02 


ph_ctg_1795 


-5.8 


4.2E-04 




Heat shock 70 kDa protein 14 


EKC26064.1 


fl_ctg_506 


4.6 


4.5E-05 


nde_ctg_5 


1.3 


4.5E-01 




Heat shock protein 68 


EKC22243.1 


fl_ctg_2557 


24.6 


1.7E-02 


ph_ctg_1 676 


15.9 


1.2E-06 




Heat shock protein 70 B2, type2 


EKC30019.1 


fl_ctg_1763 


4.61 


0 


ph_ctg_1 662 


1.66 


4.3E-21 




Heat shock protein beta-1, type 2 


EKC40046.1 


fl_ctg_683 


21.5 


1.6E-10 


ph_ctg_2231 


10.6 


1 .8E-243 




Heat shock protein beta-1, type 3 


EKC40046.1 


fl_ctg_1 192 


58.7 


5.0E-07 


ph_ctg_1 086 


21.3 


9.3E-78 




Heat shock protein HSP 90-alpha 1 


EKC25687.1 


fl_ctg_2540 


9.5 


5.5E-18 


ph_ctg_1 044 


9.5 


8.9E-65 




Stress-70 protein, mitochondrial 


EKC 18038.1 


fl_ctg_1 597 


3.4 


7.0E-03 


ph_ctg_1927 


1.6 


3.8E-08 




BAG molecular chaperone regulator 4 


EKC42633.1 


fl_ctg_2336 


27.9 


6.0 E- 16 










Calnexin 


EKC32723.1 


fl_ctg_224 


6.1 


2.4E-08 










Calreticulin 


EKC23905.1 


fl_ctg_1087 


2.7 


1.1E-16 










DnaJ-like protein subfamily B member 4 


EKC35 183.1 


fl_ctg_2051 


6.3 


1.1E-41 










Heat shock 70 kDa protein 12A, type 2 


EKC 19496.1 


fl_ctg_2338 


3.6 


2.7E-02 










Heat shock factor protein 1 


EKC33988.1 


fl_ctg_757 


6.1 


2.3E-02 










Heat shock protein 70 B2, typel 


EKC21713.1 


fl_ctg_2493 


8.9 


2.6E-03 










Heat shock protein beta-1, type 1 


EKC27576.1 


fl_ctg_1 863 


173.2 


2.9E-10 










HSPB1 -associated protein 1 


EKC31 176.1 


fl_ctg_1 565 


4.0 


1.3E-05 










Heat shock 70 kDa protein 12B 


EKC29173.1 








ph_ctg_1 045 


7.4 


3.2E-06 




Antioxidant/Oxidative Stress Response (p 


= 0.1534) 
















Calpain-5 


EKC27521.1 


fl_ctg_1914 


1.7 


1 .1 E-02 


ph_ctg_688 


4.4 


1.1E-05 




Ceruloplasmin 


EKC34134.1 


fl_ctg_810 


4.8 


6.8E-31 


ph_ctg_1 266 


-17.3 


1 .6E-02 




Dual oxidase 


EKC42615.1 


fl_ctg_941 


-7.0 


2.1E-03 


ph_ctg_2002 


-4.4 


2.6E-02 




Glutathione peroxidase 1 


EKC25678.1 


fl_ctg_2179 


2.9 


2.6E-14 


ph_ctg_1277 


2.2 


9.9E-12 




Pantetheinase 


EKC3 1982.1 


fl_ctg_1534 


2.6 


1.0E-04 


nde_ctg_ 1 1 


1.2 


6.5E-01 




Protein toll 


EKC33894.1 


fl_ctg_829 


10.9 


2.7E-04 


ph_ctg_1155 


3.1 


4.9E-02 




Selenide, water dikinase 


EKC33185.1 


nde_ctg_4 


1.3 


2.4E-02 


ph_ctg_1540 


2.3 


7.7E-20 




Selenoprotein M 


NP_001 18681 1.1 


fl_ctg_1 703 


-3.7 


3.3E-02 


ph_ctg_1 989 


-2.4 


4.2E-04 




Selenoprotein P, plasma, lb 


AAH86844.1 


fl_ctg_2541 


-2.1 


1 .3E-04 


ph_ctg_1412 


-2.1 


0 




Selenoprotein T 


XP_974477.1 


fl_ctg_718 


3.0 


6.9E-10 


ph_ctg_1 960 


2.0 


3.6E-04 




Superoxide dismutase 


EKC26875.1 


fl_ctg_1390 


-3.2 


1 .4E-04 


ph_ctg_460 


-3.4 


2.2E-02 




Thioredoxin reductase 3 


EKC27867.1 


fl_ctg_1736 


-1.8 


2.1 E-02 


ph_ctg_2174 


-2.0 


7.8E-06 




X-box-binding protein 1 


EKC27673.1 


fl_ctg_2520 


33.0 


2.0E-67 


ph_ctg_388 


2.5 


1.1E-06 




Elongation factor 1 -alpha 


EKC33063.1 


fl_ctg_1966 


2.5 


1 .4E-02 










Ubiquitin-protein ligase E3C 


EKC38417.1 


fl_ctg_1472 


9.8 


1 .5E-03 










Calcium homeostasis endoplasmic 
reticulum 


EKC18572.1 








ph_ctg_2138 


6.1 


5.9E-03 




Selenoprotein W 


NP_001 159799.1 








ph_ctg_685 


-1.7 


3.6E-07 




Cell proliferation/Apoptosis (p = 0.2148) 


AP-2 complex subunit mu-1 


EKC19414.1 


fl_ctg_415 


-2.0 


2.3E-02 


ph_ctg_784 


21.5 


1 .6E-05 




Apoptosis regulator R1 


EKC20695.1 


fl_ctg_2184 


1.7 


6.0E-03 


ph_ctg_159 


2.1 


1.7E-02 
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Table 4. Cont. 



Gene description 


NCBI Accession 
number 


P. grandis 






U. tetralasmus 






t-test 
p<0.05 






Contig ID 


FC 


p-value 


Contig ID 


FC 


p-value 




BNIP1 


EKC 19480.1 


fl_ctg_43 


2.8 


4.1 E-02 


ph_ctg_1 64 


11.0 


1.3E-13 




BNIP3 


EKC35529.1 


fl_ctg_593 


3.8 


1 .6E-07 


ph_ctg_2228 


3.2 


2.3 E- 13 




Bcl-2-like protein 1 


EKC30554.1 


fl_ctg_1908 


2.8 


1.7E-27 


ph_ctg_493 


3.4 


5.4E-06 




Cathepsin L-like cysteine proteinase 


EKC36430.1 


fl_ctg_2106 


2.3 


3.6E-08 


ph_ctg_523 


-2.1 


5.0E-02 


» 


Cold shock domain-containing protein El 


EKC28502.1 


fl_ctg_1 124 


1.8 


2.3 E- 10 


ph_ctg_2058 


2.0 


2.6E-04 




Corticosteroid 1 1-beta-dehydrogenase 2 


EKC40090.1 


fl_ctg_1 159 


32.2 


1 .2E-06 


ph_ctg_608 


7.4 


5.0E-02 


* 


Cyclin-I 


EKC18711.1 


fl_ctg_2271 


2.0 


6.2E-09 


ph_ctg_2184 


2.0 


1.7E-05 




E3 ubiquitin-protein ligase MIB2 


EKC39104.1 


fl_ctg_2263 


-6.2 


2.1 E-02 


ph_ctg_1106 


3.6 


4.6E-05 


» 


GADD45 alpha 


EKC3 1645.1 


fl_ctg_62 


28.2 


1.2E-05 


nde_ctg_ 12 


-1.2 


6.3E-01 




Inhibitor of apoptosis protein 


EKC20774.1 


fl_ctg_1444 


14.7 


7.8E-24 


ph_ctg_1 060 


11.9 


6.7E-03 




Krueppel-like factor 5 


EKC 17625.1 


fl_ctg_1 639 


25.8 


8.1E-60 


ph_ctg_1332 


3.7 


3.0E-04 




Polycomb complex protein BMI-1 


EKC23428.1 


fl_ctg_1 135 


2.7 


7.7E-06 


ph_ctg_1337 


2.9 


2.0E-08 




Protein BTG1 


EKC25710.1 


fl_ctg_1732 


9.1 


8.4E-175 


ph_ctg_167 


1.5 


1 .OE-02 




Serine/threonine-protein kinase Pim-3 


EKC 19760.1 


fl_ctg_1493 


3.7 


6.6E-13 


ph_ctg_1242 


-2.4 


1.3E-08 




Stress-induced-phosphoprotein 1 


EKC 18743.1 


fl_ctg_1646 


3.5 


3.5E-08 


ph_ctg_1651 


2.5 


9.8E-08 




TMBIM4 


EKC24718.1 


fl_ctg_853 


2.8 


1.3 E-02 


ph_ctg_2481 


2.0 


3.1E-12 




Ubiquitin 


EKC42613.1 


fl_ctg_971 


14.0 


5.0E-1 3 


ph_ctg_231 


2.0 


1.0E-31 


* 


Antigen KI-67 


EKC25595.1 


fl_ctg_1 31 6 


—4.9 


1 .3 E-02 










Caspase-10 


EKC25819.1 


fl_ctg_1 166 


8.2 


2.8E-05 










CDH1-D 


AAL31 950.1 


fl_ctg_1990 


2.4 


2.8E-06 










MAP kinase signal-integrating kinase 1 


EKC18571.1 


fl_ctg_1707 


3.9 


1.8E-10 










Protein SET 


EKC28777.1 








ph_ctg_761 


-2.6 


2.7E-04 




Putative inhibitor of apoptosis 


EKC26950.1 








ph_ctg_868 


86.8 


6.5E-26 




Ubiquitin-conjugating E2 variant 1 
(UEV1A) 


EKC31486.1 








ph_ctg_910 


4.5 


5.1E-14 




Energy Metabolism (p = 0.1482) 


5 -AMP-activated protein kinase beta-2 


EKC39833.1 


fl_ctg_2083 


-4.3 


1 .2E-05 


ph_ctg_2306 


-2.1 


2.9E-02 




ATP synthase F0 subunit 6 


YP_003355006.1 


fl_ctg_2323 


-7.1 


5.8E-08 


nde_ctg_6 


3.3 


2.7E-01 




Cholecystokinin receptor 


EKC30410.1 


fl_ctg_1030 


3.8 


1 .1 E-04 


ph_ctg_337 


-3.3 


4.9E-02 




Cytochrome b 


YP_003354999.1 


fl_ctg_1821 


-5.6 


0 


ph_ctg_589 


-2.7 


2.1 E-02 




Cytochrome c oxidase subunit 1 


YP_003354995.1 


fl_ctg_1261 


-3.4 


0 


ph_ctg_288 


-1.6 


1.3E-13 


* 


Cytochrome c oxidase subunit II 


YP_003354996.1 


fl_ctg_1 138 


-3.3 


0 


ph_ctg_23 


-1.9 


3.1 E-04 




Cytochrome c oxidase subunit III 


YP_003355007.1 


fl_ctg_1487 


-4.5 


0 


nde_ctg_7 


-2.4 


3.1E-01 




Glutamine synthetase 2 cytoplasmic 


EKC27195.1 


fl_ctg_1817 


5.3 


1 .5E-24 


ph_ctg_1 688 


2.0 


1.6E-24 




Group XIIA secretory phospholipase A2 


EKC23468.1 


fl_ctg_2457 


2.8 


1 .6E-03 


ph_ctg_1 1 98 


10.9 


5.0E-05 




NADH dehydrogenase subunit 1 


YP_003355001.1 


fl_ctg_1537 


-6.0 


0 


nde_ctg_ 14 


-1.3 


8.8E-01 




NADH dehydrogenase subunit 4 


YP_003355003.1 


fl_ctg_1864 


-2.0 


1 .5E-03 


nde_ctg_8 


-1.1 


9.5E-01 




NADH dehydrogenase subunit 5 


YP_003355000.1 


fl_ctg_1887 


-3.8 


2.5E-12 


ph_ctg_503 


1.7 


1.5E-03 




Succinyl-CoA ligase 


EKC28914.1 


fl_ctg_1283 


-2.2 


2.7E-03 


ph_ctg_1091 


-3.5 


6.2E-06 




ATF5 activating transcription factor 5 


EKC37224.1 


fl_ctg_360 


2.4 


9.5E-07 










Immune (p = 0.7569) 


B-cell lymphoma 3-encoded protein 


EKC39819.1 


fl_ctg_200 


6.2 


1.2E-02 


ph_ctg_1522 


-2.4 


4.0E-06 




Histone H4 transcription factor 


EKC36777.1 


fl_ctg_2257 


6.0 


5.8E-08 


nde_ctg_9 


1.9 


5.5E-01 




Interferon regulatory factor 2 


EKC43 156.1 


fl_ctg_1 121 


-2.5 


1 .4E-06 


ph_ctg_76 


-4.2 


3.3E-10 




lnterleukin-1 receptor-associated 
kinase 1 BP1 


AER35157.1 


nde_ctg_2 


1.1 


6.0E-01 


ph_ctg_1 648 


-1.7 


1.3 E-02 




Lipopolysaccharide-binding protein 


EKC33909.1 


fl_ctg_2167 


-2.8 


9.7E-03 


ph_ctg_1475 


-2.0 


1.2E-03 




LPS-induced tumor necrosis 
factor-alpha 


EKC36231.1 


fl_ctg_1784 


-4.3 


2.0E-03 


ph_ctg_1823 


-2.3 


4.9E-02 
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Table 4. Cont. 



Gene description 


NCBI Accession 
number 


P. grandis 






U. tetralasmus 






t-test 
p<0.05 






Contig ID 


FC 


p-value 


Contig ID 


FC 


p-value 




Lysozyme 


EKC39289.1 


fl_ctg_477 


2.1 


4.1E-02 


ph_ctg_2014 


2.3 


3.0E-10 




NLR family, pyrin domain containing 3 


XPJD03730107.1 


fl_ctg_758 


4.5 


9.8E-04 


ph_ctg_1933 


2.7 


1.5E-03 




NF-kappa-B inhibitor cactus 


EKC37718.1 


fl_ctg_2206 


3.9 


1.9E-18 


ph_ctg_935 


2.1 


2.5E-04 


» 


Peptidoglycan-recognition protein SC2 


EKC42861.1 


fl_ctg_1392 


3.4 


3.0E-03 


ph_ctg_1 1 82 


3.4 


1.6E-13 




Phospholipid scramblase 2, partial 


EKC25020.1 


fl_ctg_585 


-4.0 


2.0E-02 


ph_ctg_894 


14.3 


5.4E-06 




Serine/threonine-protein kinase TBK1 


EKC41 453.1 


fl_ctg_1373 


7.0 


3.2E-31 


ph_ctg_1 1 1 5 


11.3 


3.4E-03 




Slit-like protein 1 protein 


EKC33230.1 


fl_ctg_1 877 


6.5 


7.0E-06 


nde_ctg_ 13 


1.2 


5.2E-01 




TNFAIP3-interacting protein 2 


EKC38434.1 


fl_ctg_840 


4.7 


4.3E-04 


ph_ctg_1 846 


2.5 


2.2E-02 




TNF receptor-associated factor 2 


bKLJlboz.l 


nde_ctg_ 7 


1.9 


5.4E-0 7 


ph_ctg_783 


4.5 


1.7E-03 




Toll-like receptor 2 type-2 


bKLloJJ/.l 


fl_ctg_1 332 


3.8 


9.9E-06 


ph_ctg_1 237 


4.0 


1 .3E-1 0 




Toll-like receptor 3 


bKLJbybo.l 


Tl_ctg_l 369 


—3.0 


3.4E-03 


ph_ctg_26 


4.6 


7.1 E-03 




14-3-3 protein epsilon 


EKL29146.1 


Tl_ctg_1686 


—4.2 


9.8E-1 3 










Myeloid differentiation primary 
response 88 


bKL40UOD.l 


TI_ctg_23/8 


3.6 


3.0E-03 










Taxi -binding protein 1 -like protein B 


EKC20118.1 


fl_ctg_2123 


-1.5 


8.1E-03 










14-3-3 protein zeta 


EKC18419.1 








ph_ctg_1953 


-7.8 


1.3E-04 




Complement C3 


EKC24393.1 








ph_ctg_2299 


-1.6 


1.7E-06 




Defensin 


ADO32580.1 








ph_ctg_2420 


100.1 


1.9E-06 




EIF4EBP1 


EKC38736.1 








ph_ctg_889 


2.5 


1.1E-07 




NLR family, pyrin domain containing 1 


XPJD02731351.1 








ph_ctg_1301 


13.2 


2.2E-03 




TNF receptor-associated factor 3 


EKC22057.1 








ph_ctg_2244 


3.3 


4.5 E-03 




Toll-like receptor 13 


EKC26527.1 








ph_ctg_1671 


-10.7 


1 .8E-02 




Cytoskeletal (p = 0.01 49) 


Collagen alpha-2(l) chain 


EKC 19703.1 


fl_ctg_1 852 


-2.3 


1 .OE-02 


ph_ctg_1792 


-1.9 


1.5E-05 




Fibrillin! -like 


EKC40135.1 


fl_ctg_1999 


-2.3 


2.0E-02 


ph_ctg_2292 


-6.5 


1.2E-03 




Microtubule-associated protein 2 


EKC35 104.1 


fl_ctg_1783 


40.1 


4.4E-08 


ph_ctg_664 


5.9 


6.1E-04 




Myosin catalytic light chain LC-1 


EKC20656.1 


nde_ctg_3 


1.2 


8.4E-02 


ph_ctg_1 040 


-13.0 


2.1E-04 




Protocadherin Fat 4 


EKC27752.1 


fl_ctg_2135 


1.9 


8.5E-03 


ph_ctg_1 385 


-16.5 


2.8E-07 




Ras-like GTP-binding protein Rhol 


EKC37733.1 


fl_ctg_982 


15.5 


3.8E-03 


ph_ctg_997 


-2.7 


4.9E-02 




Tubulin alpha-IC chain 


EKC 19808.1 


fl_ctg_738 


-2.4 


1.3E-10 


ph_ctg_1023 


-2.2 


5.4E-09 




Tubulin beta chain 


EKC30036.1 


fl_ctg_408 


33.1 


7.9E-31 


ph_ctg_2439 


3.6 


1.2E-03 


* 


Dynein heavy chain 6, axonemal 


EKC31 674.1 


fl_ctg_2151 


-2.8 


7.6E-04 










Tubulin alpha chain 


EKC22339.1 


fl_ctg_442 


-3.4 


3.7E-08 










Other (p = 0.8545) 


Chitin synthase 1 


EKC23809.1 


fl_ctg_835 


-3.5 


2.1E-04 


ph_ctg_1713 


-4.3 


3.5E-06 




Chloride intracellular channel exc-4 


EKC30516.1 


fl_ctg_1892 


5.2 


6.5E-22 


nde_ctg_10 


1.4 


1.8E-02 




Far upstream element-binding 
protein 3 


NP_001 034426.1 


fl_ctg_614 


-31.2 


0 


ph_ctg_666 


-1.9 


1.8E-02 




Fibulin-2 


EKC31 789.1 


fl_ctg_604 


-3.9 


6.4E-03 


ph_ctg_1457 


-3.8 


5.0E-08 




Hemicentin-1 


EKC33127.1 


fl_ctg_1013 


15.8 


4.9E-21 


ph_ctg_1370 


-2.4 


0 




Krueppel-like factor 1 1 


EKC40859.1 


fl_ctg_1 120 


3.9 


4.1E-07 


ph_ctg_245 


3.6 


5.0E-02 




Outer dense fiber protein 3 


NP_001 171751.1 


fl_ctg_846 


-12.2 


8.1E-04 


ph_ctg_803 


4.0 


4.9E-02 




Putative accessory gland protein 


ABG01813.1 


fl_ctg_498 


21.0 


5.2E-20 


ph_ctg_280 


2.2 


4.5E-05 




Tribbles-like protein 2 


EKC22747.1 


fl_ctg_1325 


10.4 


2.3E-164 


ph_ctg_2106 


1.9 


3.0E-09 




Water and ammonia transporting 
aquaporin 


EKC35418.1 


fl_ctg_1403 


2.3 


1 .5E-04 


ph_ctg_1774 


-1.6 


2.4E-05 




Epithelial membrane protein 2 


EKC25041.1 


fl_ctg_113 


17.2 


3.2E-07 










Liver stage antigen 3 precursor 


XP_002257838.1 


fl_ctg_642 


-19.1 


1.3E-15 










Semenogelin II precursor 


AB052944.1 


fl_ctg_628 


-16.5 


0 
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Table 4. Cont. 





Gene description 


NCBI Accession 
number 


P. grandis 






U. tetralasmus 




t-test 
p<0.05 






Contig ID 


FC 


p-value 


Contig ID 


FC p-value 




Dexrasl 


XP_002430253.1 








ph_ctg_1269 


24.8 1.7E-02 





Genes were organized in broad functional categories. Absence of contig, fold change, and p-value information for a given gene indicates that no orthologous sequence 
was detected in that particular species. Italicized contig identifiers indicate that the gene did not fit criteria for significance (^1.5 fold change, p<0.05) in the particular 
species. * indicates significantly different fold changes (FC) between species for a given gene based on a t-test of replicated fold change values (p<0.05). Differences in 
fold changes between species in each category were assessed for significance using a paired t-test (p<0.05) with values given next to the category titles. Additional 
information for each contig is available in Table S3. 
doi:1 0.1 371 /joumal.pone.0089481 .t004 



both species and a gene found to be up-regulated uniquely in each 
species (Figure 2). We found that heat shock gene responses in P. 
grandis showed faster and higher levels of induction than observed 
in U. tetralasmus. All four shared heat shock proteins (CRYAB, 
HSP90AA1, HSP70B2-2, HSPB1-3) showed significant induction 
of expression in the more sensitive P. grandis by 24 h, while up- 
regulation of these genes in U. tetralasmus at the same timepoint was 
either relatively modest or not observed. Expression of the four 
shared genes trended higher throughout the heat and desiccation 
timepoints, with the exception of HSPB1-3 in P. grandis which 
peaked at 24 h and then declined. Similar patterns of up- 
regulation were also observed for BAG molecular chaperone 
regulator 4 (BAG4) in P. grandis and defensin (DEF) in U. 
tetralasmus. All examined genes had declined from their 72 h 
expression levels in the recovery timepoint, with most genes 
approaching homeostatic (0 h) levels. Fold changes (relative to 0 h) 
were significantly higher in P. grandis than U. tetralasmus at most 
tested timepoints (Figure 2). Based on this data subset, it seems 
likely that that other key genes identified at 72 h also responded 
sensitively to heat/desiccation at earlier timepoints and returned 
to normal levels as experimental stressors were relaxed (recovery 
period). 

Discussion 

Recent studies of mussel behavioral and physiological responses 
to drought [2,3] have revealed an array of strategies to maximize 
survival, growth, and reproduction. Mussel species inhabiting 
small streams during drought periods can be stranded by receding 
waters, forcing them to burrow into cooler, wetter sediments or to 
track waters to deeper pools, when available. These differing 
behavioral approaches (tracking, burrowing, or tracking then 
burrowing) likely result in differential fluctuations in population 
sizes and stream bed distributions among species based on annual 
variations in the hydrologic cycle. Physiologically, stranded 
mussels face the simultaneous stressors of increased temperature 
and tissue desiccation which can, ultimately, impact survival but, 
more immediately, regulate the evolved avoidance strategies of 
each species. We sought here to examine the transcriptional 
regulation of differing tolerances to drought in two mussel species 
representing physiological and behavioral extremes. U. tetralasmus 
(pondhorn) can live up to 2 years out of water at cool 
temperatures, while P. grandis (floater) can tolerate only a few 
weeks of desiccation [30]. Recent field and laboratory studies by 
our group revealed burrowing behavior in U. tetralasmus and 
tracking behavior in P. grandis [3]. In this study, we carried out 
RNA-Seq transcriptome profiling on P. grandis and U. tetralasmus 
exposed to experimental conditions simulating those of a drought- 
induced stranding event, i.e. heat and desiccation. 



Prior to initiation of this research, few genetic resources were 
available for either P. grandis or U. tetralasmus, with only a handful 
of mitochondrial sequences currently in public databases. There- 
fore, prior to examining transcriptional responses in either species, 
we needed to sequence and assemble high-quality transcriptomes. 
Use of Trinity [18] provided superior assembly when compared 
with Trans-ABySS and generated contigs with average lengths 
greater than 800 bp. Although it is difficult to predict gene 
numbers based on comparison with other invertebrate species 
known for high rates of mutation and diversification [31], 
estimates of gene numbers from the draft genomes of the Pacific 
oyster (Crassostrea gigas) and the pearl oyster (Pinctada Jiicata) are 
available. These bivalve species encode 28,027 and 23,257 genes, 
respectively [29,32], compared with 22,616 and 20,610 predicted 
unigene annotations in P. grandis and U. tetralasmus, respectively 
(Table 2). Based on predicted low rates of sequence conservation 
across species, we estimate that the P. grandis and U. tetralasmus de 
novo assemblies captured at least 70% of the coding genes of each 
species. 

Sampled tissues and pooling strategies inevitably impact the 
size, scope, and nature of captured transcriptomes and associated 
downstream expression profiles. Here, we sampled tissues from the 
foot of individual mussels using a biopsy punch and pooled 
individual samples to construct replicate pools for RNA-Seq 
analysis as in previous studies [15,26-28]. While pooling may 
mask some of the individual variation in homeostatic and induced 
gene expression profiles, it also serves to normalize outliers and 
direct focus to key shared responses. Putative biomarkers 
developed here can be field-validated on individual samples in 
future studies. 

The mussel foot is a tough, extendible muscular organ used 
primarily for locomotion and, in juveniles, for anchoring the 
animal to substrate via byssal threads. No reports of gene 
expression in the foot of unionid mussels are available to-date, 
although flow regimen has recently been reported to impact byssal 
thread production in juvenile unionids [33]. In the zebra mussel, 
Dreissena poljmorpha, gene and protein expression of the foot organ 
have been studied in the context of adhesion [34—36]. Our 
primary consideration in choosing to examine the foot transcrip- 
tome was to obtain samples rapidly and in a non-lethal manner 
[37]. Given that 70% of freshwater mussels are currently imperiled 
[4], future wide-scale assessments of individual or population 
health in the field will require sampling techniques that are rapid 
and minimally invasive. We, therefore, wanted to obtain relevant 
signatures of heat/ desiccation stress from a tissue available without 
sacrifice or extensive manipulation of the organism. Additionally, 
expression signatures in the foot may reflect the resulting 
differential behavioral impulses for locomotion, adhesion and/or 
burrowing. While beyond the scope of the present study, a 
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Figure 1 . Comparison of fold changes between RNA-Seq and QPCR results in selected genes of P. grandis (A) and U. tetralasmus (B) at 
72 h of heat/desiccation exposure. QPCR fold changes are relative to control samples and normalized by changes in beta-actin values. Gene 
abbreviations are: Defensin, DEF; Inhibitor of apoptosis protein, IAP; Heat shock protein beta-1, type 3,HSPB1-3; Alpha-crystallin B chain, CRYAB; NLR 
family, pyrin domain containing 1, NLRP1; Heat shock protein HSP 90-alpha 1, HSP90AA1; heat shock 70 kDa protein 12B-like, HSPA12B; BCL2/ 
adenovirus E1B interacting protein 1 -like, BNIP1; Heat shock protein 70 B2, type2, HSP70B2-2; Kruppel-like factor5, KLF5; Calpain 5, CAPN5; Toll-like 
receptor 13, TLR13; Heat shock protein beta-1, type 1,HSPB1-1; BAG family molecular chaperone regulator 4, BAG4; Heat shock protein 70 B2,type1, 
HSP70B2-1; Liver stage antigen 3 precursor, LSA-3. Contig IDs are available in Table S1. 
doi:1 0.1 371 /journal.pone.0089481 .g001 



comprehensive tissue expression atlas of certain representative, 
non-endangered mussel species should be considered in the future. 

Several transcriptomic and proteomic studies in bivalves over 
the last decade have expanded our understanding of the 
components of thermal stress responses in these species. The most 
comprehensive studies have been carried out in marine mussels 
[10-13] and in the Pacific oyster [8,9,29]. With the exception of 
Zhang et al. 2012 [30], other studies have utilized microarray 
technology, potentially missing rare, previously un-sequenced 
transcripts. Recently, we conducted a smaller-scale heat stress 
study in the unionid mussel Villosa lienosa utilizing RNA-Seq [15]. 
Examined together, several key components of bivalve heat stress 



responses can be seen in these studies. These include changes in 
molecular chaperones, oxidative stress responses, changes in cell 
turnover and death, changes in energy metabolism, immune and 
inflammatory responses, and cytoskeletal reorganization. We 
observed similar dysregulation of genes involved in these pathways 
in our datasets from P. grandis and U. tetralasmus (Table 4). Below 
we highlight putative functions of selected key genes from these 
categories which may underlie shared and species-specific 
responses to drought. 

Chaperones/Heat Shock Proteins. Up-regulation of mo- 
lecular chaperones, particularly HSP70 family members, has long 
been understood as a classic indicator of protein damage due to 
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Figure 2. Temporal profiling of key genes by QPCR in P. grandis and U. tetralasmus. Fold changes are relative to control values (set as 1) and 
normalized by changes in beta-actin values. * indicates significantly different fold changes between species (p<0.05). Error bars reflect standard error 
of the mean. Gene abbreviations are: Alpha-crystallin B chain, CRYAB; Heat shock protein HSP 90-alpha 1, HSP90AA1; Heat shock protein 70 B2, type2, 
HSP70B2-2; Heat shock protein beta-1, type 3,HSPB1-3; BAG family molecular chaperone regulator 4, BAG4; Defensin, DEF. Contig IDs are available in 
Table SI. 

doi:1 0.1 371 /journal.pone.0089481 .g002 



heat [12,38]. Heat shock proteins (HSPs) interact with stress- 
denatured proteins, preventing their aggregation. An array of 
molecular chaperones was induced in both species by exposure to 
72 h of elevated temperatures and desiccation (Table 4). Only one 
chaperone (HSP70-12A1) was down-regulated in both species. A 
general pattern of higher expression of shared chaperones in the 
less-tolerant mussel species was present, with 8 of 15 shared 
chaperones showing a significandy higher fold increase by P. 
grandis compared to U. tetralasmus, but none showing a significantly 
lower fold increase. For example, a small heat shock protein, 
HSPBl-type 3 was induced 58.7-fold in P. grandis compared with 
2 1 .3-fold in U. tetralasmus. An additional nine HSP family members 



were unique to P. grandis, including HSPBl-type 1 (173.2-fold). 
These differences were confirmed in selected HSP genes analyzed 
by QPCR (Figure 1). Further temporal profiling of four shared 
molecular chaperones showed that species differences in chaper- 
one expression were not confined to the 72 h timepoint but were 
also present at 24 h and 48 h. Uniomerus tetralasmus HSPs showed 
limited induction prior to 72 h, while P. grandis HSPs were rapidly 
induced by 24 h. Induction of the HSP response was temporal in 
nature, however, as expression levels of the tested genes all had 
returned to baseline levels by 5 d after cessation of the emersion/ 
heat exposure (Figure 2). In previous comparisons of marine 
mussel species and populations with differing heat tolerances, it 
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was observed that differences in small HSP (HSPB-type) up- 
regulation were among the most significant signatures of heat 
tolerance/sensitivity [10-13]. There, small HSPs were induced at 
lower temperatures and/ or to higher levels in heat-sensitive mussel 
groups. Small HSPs, often denoted as HSPBs, are a family of heat 
shock proteins involved in regulation of apoptosis, oxidative stress, 
and cytoskeletal regulation [39,40]. Here also, HSPB family 
members were highly induced in P. grandis relative to U. tetralasmus 
indicating their potential importance in determining drought 
tolerance in unionid mussels. 

Cell Proliferation/Apoptosis. Accumulation of oxidative 
damage can induce a number of pro-apoptotic signaling pathways. 
Conversely, protective/adaptive mechanisms can seek to block 
excessive apoptotic signals and stimulate cell proliferation and a 
return to homeostasis [41,42]. Many of the components of these 
processes were observed in the sets of differentially expressed genes 
from both mussel species (Table 4). However, differential 
regulation of several genes may indicate crucial differences 
impacting drought tolerance. Genes needed for protein degrada- 
tion such as ubiquitin were higher in the sensitive P. grandis (14.0- 
fold up-regulated) than in U. tetralasmus (2.0-fold). GADD45, a key 
indicator of DNA damage, was up-regulated 28.2-fold in P. grandis 
but was not induced in U. tetralasmus. GADD45 was similarly up- 
regulated in Mytilus galloprovincialis following exposure to metal salts 
[43]. Similarly, caspase-10, a component of the execution phase of 
cell apoptosis was up-regulated 8.2-fold in P. grandis but unchanged 
in U. tetralasmus. Interestingly, an inhibitor of apoptosis (LAP) gene 
showed the second highest up-regulation (86.8-fold) of all 
differentially expressed genes in the drought-tolerant U. tetralasmus. 
No clear IAP orthologue was identified from P. grandis. Lastly, 
ubiquitin-conjugating enzyme E2 variant 1 (UEV1A) was differ- 
entially expressed only in U. tetralasmus (Table 4). In mammals, IAP 
genes interact with UEV1A and other ubiquitin ligases such as 
MIB2 (also differentially expressed here) to inhibit stress-induced 
apoptosis and facilitate cell proliferation/survival [44-46]. These 
pathways are mediated through NF-kB activation and, therefore, 
often dovetail with modulation of components of the innate 
immune system [47,48_ENREF_50] . Taken together, the sensitive 
P. grandis showed induced pro-apoptotic signaling pathways and 
indicators of DNA damage, whereas the tolerant species U. 
tetralasmus showed greater inhibition of apoptosis and little 
indication of DNA damage. These molecular differences are likely 
an important component of the observed physiological differences 
in drought tolerance between the two species. Further studies are 
needed in molluscan species to confirm the conserved function of 
these components in regulation of cell survival. 

Immune. One area of concern for drought-exposed mussels 
is that stress responses induced by emersion and/or rising air/ 
water temperatures will impact immune health increasing 
susceptibilities to pathogens encountered in their environment 
[9,33]. Indeed, we observed the dysregulation of a number of 
important mediators of immunity in both mussel species (Table 4). 
These included the down-regulation of innate cytokines such as 
TNF-a and the down-regulation of pathogen-recognition recep- 
tors such as Toll-like receptor 3. Several components of the NF-kB 
signaling pathway including TRAF2, TRAF3, and TBK1 were 
either present only in U. tetralasmus or were induced to a greater 
extent, potentially regulating cell survival mechanisms as discussed 
above. 

Of particular interest in the immune category was the strong up- 
regulation of a U. tetralasmus defensin (up 100.1 -fold; highest FC 
among pondhorn genes). An orthologous defensin was not 
identified among differentially expressed genes in P. grandis. 
Defensins are a family of peptides with antimicrobial activities 



against a range of bacterial and fungal pathogens. While their 
traditional immune roles have been characterized in marine 
mussels and oysters [49_ENREF_51,50], little is known about 
their functions in freshwater mussel. Recendy, however, Xu and 
Faisal [51] described a defensin from zebra mussel, D. polymorpha, 
with antimicrobial activity. Notably, zebra mussel defensin 
expression was the highest in the foot of all tested tissues and 
cell types and expression levels corresponded with status of 
byssogenesis. The authors there speculated that defensin may be 
providing biological protection against microbial degradation of 
byssal threads or may be playing a more fundamental role in 
byssogenesis. Here, in adult unionids, observed differences in 
defensin up-regulation may reflect differential use of the foot in 
burrowing. The burrowing species U. tetralasmus may scratch or 
injure its foot in burrowing, resulting in protective secretion of the 
antimicrobial defensin. By QPCR analysis, up-regulation of U. 
tetralasmus defensin was even higher than that indicated by RNA- 
Seq (Figure 1) and expression levels rose with the duration of the 
challenge (Figure 2). Future studies should examine potential 
functions of unionid defensins in the context of immune status, 
burrowing, byssal thread production, and drought tolerance. 

Cytoskeletal. The buildup of free radicals associated with 
cellular stress leads to damage of the actin cytoskeleton. Changes 
to actin filament dynamics impact junctional integrity, cell cycling, 
and cellular movement. Gene and protein changes in cytoskeletal 
elements have been reported previously in response to heat stress 
in marine mussels [11,12]. In response to heat stress, cells 
modulate expression of these elements and enhance production of 
small HSPs to combat these effects [52]. In addition to the higher 
expression of small HSPs in P. grandis, larger fold changes were 
seen in genes associated with assembly of microtubules, one of the 
fibers composing the cytoskeleton. For example, tubulin beta chain 
was up-regulated 33.1 -fold in P. grandis, but only 3.6-fold in U. 
tetralasmus (Table 4). These differences are another indication of 
higher levels of cellular stress in P. grandis. 

Other. We noted a handful of additional genes with putative 
functions outside one of the larger response categories. These 
included genes with functions in ion transport, reproduction, shell 
structure, and circadian rhythms (Table 4). Interestingly, in both 
species drought stress led to reduced expression of chitin synthases, 
key enzymes involved in synthesis of bivalve shells, indicating 
potential impacts on growth and health from prolonged drought 
exposure even in "tolerant" species [53]. Finally, dexrasl, a gene 
which, in mammals, is responsible for transducing signals 
maintaining circadian rhythms of behavior and physiology [54] , 
was sharply induced in U. tetralasmus alone (24.8-fold up- 
regulation). Given the ability of stranded mussels of this species 
to endure long periods of emersion, it would be of great interest to 
determine whether dexrasl is responsible for the biological 
reprogramming needed for tolerating harsh environmental condi- 
tions. 

Conclusions 

Freshwater mussels in the southeastern US are imperiled by 
increasingly frequent drought conditions. Here we examined the 
molecular underpinnings of differing behavioral and physiological 
responses to drought in two co-occurring mussel species. RNA-Seq 
transcriptome profiling successfully captured and assembled high 
quality transcriptomes for P. grandis and U. tetralasmus and 
facilitated comparison of gene expression profiles following 
heat/desiccation exposure. Shared molecular responses included 
changes in molecular chaperones, oxidative stress profiles, 
cell cycling, energy metabolism, immunity, and cytoskeletal 
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rearrangements. Significantly higher induction of molecular 
chaperones and cytoskeletal elements indicated higher levels of 
cellular stress in P. grandis while our results also highlighted 
individual genes potentially contributing to drought tolerance in 
U. tetralasmus. From this foundation, future studies should seek to 
validate putative biomarkers revealed here in individual mussels 
from P. grandis and U. tetralasmus collected in field settings, as well 
as examining their utility in assessing stress levels in a broader 
array of unionid species. 
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Figure SI Gene ontology (GO) term categorization and 
distribution of assembled Trinity contigs encoding genes in P. 
grandis (A) and U. tetralasmus (B). 
(TIF) 
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(DOCX) 

Table S2 Summary of Illumina expressed short reads produc- 
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Paired-end reads were generated on a HiSeq 2000 instrument. 
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